This notebook aims to independently calculate the dose based on the dwell positions within an Oncentra® RT Dicom Plan, and compare the resulting dose grid to an exported Oncentra® dose grid within an RT Dicom Dose file.
Collect a range of brachytherapy dicom files that are able to be placed within the Gihub repository that can be used for testing. Aim to support as many brachy dicom formats as possible.
Once egs_brachy (https://doi.org/10.1088/0031-9155/61/23/8214) is available for use I would like to directly implement that model based dose calculation method to be able to provide checking for both TG43 and Monte Carlo based algorithms.
It is also planned that upper and lower 95% confidence interval doses will be reported when uncertainties due to catheter movement, catheter reconstruction, and calculation uncertainties are taken into account. Using the structure dicom file these 95% confidence interval doses can be converted to comparative DVHs.
The final stage is for this code to rerun a dwell time optimisation to see if any plan improvement is possible and compare the robustness of the original plan to positioning uncertainties with the calculated comparative plan.
Copyright (C) 2016-2017 Simon Biggs
This program is free software: you can redistribute it and/or modify it under the terms of the GNU Affero General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more details.
You should have received a copy of the GNU Affero General Public License along with this program. If not, see http://www.gnu.org/licenses/.
Brachytherapy dicom files make use of a large number of private dicom tags. What some of these tags mean needs to be reverse engineered at times. This program has as of yet only been tested with Oncentra® dicom files, not BrachyVision™ dicom files. It has only been tested with a small subset of Oncentra® dicom files, with some of these files for certain configurations this code does not yet work. Be sure when using this code to investigate the testing figures produced to confirm that they represent what is expected within the plan.
Pay particular attention to x,y,z definitions, catheter definitions, and source orientation.
In [1]:
import numpy as np
import pandas as pd
import xarray as xr
from scipy.interpolate import RegularGridInterpolator
from copy import copy
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.collections import LineCollection
import matplotlib.gridspec as gridspec
%matplotlib inline
import itertools
from IPython.display import display, Markdown
In [2]:
source_date_directory = 'source_data/estro_website/'
# I couldn't seem to get Granero's data to work well
# source_date_directory = 'source_data/granero_paper/'
# I would like to also try Taylor's data
In [3]:
dose_rate_constant = np.squeeze(
pd.DataFrame.from_csv(source_date_directory + 'Lambda.csv', index_col=None, header=None).values)
dose_rate_constant
Out[3]:
In [4]:
# When using the Estro data it appears that a length of 0.36 agrees better than a length of 0.35.
# I am worried that during the consensus data merging that data lookups that were based on a length of 0.36
# were potentially mixed with data lookups based on a length of 0.35.
length = np.squeeze(
pd.DataFrame.from_csv(source_date_directory + 'L.csv', index_col=None, header=None).values)
# length = 0.36
length
Out[4]:
In [5]:
radial_function_table = pd.DataFrame.from_csv(source_date_directory + 'g_L(r).csv')
radial_function_data = xr.DataArray(
np.squeeze(radial_function_table.values),
coords=[
radial_function_table.index.values.astype(float)],
dims=['radius_cm'])
radial_function_data
Out[5]:
In [6]:
anisotropy_function_table = pd.DataFrame.from_csv(source_date_directory + 'F(r,theta).csv')
anisotropy_function_data = xr.DataArray(
anisotropy_function_table.values,
coords=[
anisotropy_function_table.index.values.astype(float),
anisotropy_function_table.columns.values.astype(float)],
dims=['theta_deg', 'radius_cm'])
anisotropy_function_data
Out[6]:
In [7]:
QA_along_away_table = pd.DataFrame.from_csv(source_date_directory + 'QA_along_away.csv')
QA_along_away_data = xr.DataArray(
QA_along_away_table.values.astype(float),
coords=[
QA_along_away_table.index.values.astype(float),
QA_along_away_table.columns.values.astype(float)],
dims=['z_cm', 'y_cm'])
QA_along_away_data
Out[7]:
TG43U1S1 states:
...a log-linear function for $g(r)$ interpolation was recommended in the 2004 AAPM TG-43U1 report. ... The log-linear interpolation should be performed using data points immediately adjacent to the radius of interest.
...
Appendix C 3 of the 2004 AAPM TG-43U1 report clearly specified an extrapolation method (nearest neighbor or zeroth order) for 1D dose rate distributions when $r<r_{min}$ for $g(r_{min})$. Due to the great variability in $g(r)$ based on choice of L and features of source construction, use of nearestneighbor or zeroth-order data is still recommended for extrapolation of $g_L(r)$ for $r<r_{min}$.
...
The AAPM now recommends adoption of a single exponential function based on fitting $g_L(r)$ data points for the largest two consensus $r$ values ... Specifically, a log-linear extrapolation ...
In [8]:
plt.scatter(
radial_function_data.radius_cm,
radial_function_data)
plt.xlabel("Radius (cm)")
plt.ylabel("Radial dose function")
Out[8]:
In [9]:
plt.scatter(
np.log(radial_function_data.radius_cm[1::]),
radial_function_data[1::])
plt.xlabel("$ln$( Radius (cm) )")
plt.ylabel("Radial dose function")
Out[9]:
In [10]:
radial_function_interpolation = RegularGridInterpolator(
(np.log(radial_function_data.radius_cm),),
radial_function_data,
bounds_error=False,
fill_value=None
)
def radial_function(radius):
interpolation = radial_function_interpolation(np.log(radius))
too_close_ref = radius < np.array(radial_function_data.radius_cm[0])
interpolation[too_close_ref] = radial_function_data[0]
return interpolation
r_test = np.linspace(0.03, 30, 10000)
y_test = radial_function(r_test)
plt.figure()
plt.scatter(
np.log(radial_function_data.radius_cm),
radial_function_data)
plt.plot(
np.log(r_test),
y_test)
plt.xlabel("$ln$( Radius (cm) )")
plt.ylabel("Radial dose function")
plt.figure()
plt.scatter(
radial_function_data.radius_cm,
radial_function_data)
plt.plot(
r_test,
y_test)
plt.xlabel("Radius (cm)")
plt.ylabel("Radial dose function")
Out[10]:
TG43U1S1 states:
"...linear-linear interpolation method for $F(r,\theta)$, as a function of $r$ and $\theta$ is appropriate, and should be based on the two data points for each variable located immediately adjacent to the interpolated point of interest."
"...the nearestneighbor or zeroth-order approach presented in Appendix C of the 2004 AAPM TG-43U1 report is still recommended for $F(r,\theta)$ extrapolation for $r < r_{min}$ and also for $r > r_{max}$."
In [11]:
radius_cm_mesh, theta_deg_mesh = np.meshgrid(
anisotropy_function_data.radius_cm,
anisotropy_function_data.theta_deg)
plt.scatter(
radius_cm_mesh,
theta_deg_mesh,
c=anisotropy_function_data,
s=2)
c = plt.colorbar()
plt.xlabel("Radius (cm)")
plt.ylabel("Theta (degrees)")
c.set_label("Anisotropy function")
In [12]:
find_nan = np.where(np.isnan(anisotropy_function_data))
row = find_nan[0]
column = find_nan[1]
unique_row = np.unique(row)
for row_val in unique_row:
ref = row == row_val
last_nan = np.max(column[ref])
given_row = anisotropy_function_data[row_val, :]
nan_ref = np.isnan(given_row)
anisotropy_function_data[row_val, nan_ref] = given_row[last_nan + 1]
In [13]:
plt.scatter(
radius_cm_mesh,
theta_deg_mesh,
c=anisotropy_function_data,
s=2)
c = plt.colorbar()
plt.xlabel("Radius (cm)")
plt.ylabel("Theta (degrees)")
c.set_label("Anisotropy function")
In [ ]:
In [14]:
marker_cycle = itertools.cycle((',', 'o', '*', '^'))
colour_list = [
item['color']
for item in mpl.rcParams['axes.prop_cycle']
]
colour_cycle = itertools.cycle(colour_list)
for data in anisotropy_function_data.T:
plt.scatter(
data.theta_deg, data, label="radius = {} cm".format(np.array(data.radius_cm)),
marker=next(marker_cycle), c=next(colour_cycle))
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Theta (degrees)')
plt.ylabel('Anisotropy function')
Out[14]:
In [15]:
anisotropy_theta_deg = np.array(anisotropy_function_data.theta_deg)
anisotropy_radius_cm = np.array(anisotropy_function_data.radius_cm)
anisotropy_function_data_array = np.array(anisotropy_function_data)
anisotropy_function_linear_interpolation = RegularGridInterpolator(
(anisotropy_theta_deg, anisotropy_radius_cm), anisotropy_function_data_array,
bounds_error=False,
fill_value=np.nan)
anisotropy_function_nearest_interpolation = RegularGridInterpolator(
(anisotropy_theta_deg, anisotropy_radius_cm), anisotropy_function_data_array,
bounds_error=False,
fill_value=None,
method='nearest')
def anisotropy_function(radius, theta):
points = np.vstack([theta, radius]).T
linear_interpolation = anisotropy_function_linear_interpolation(points)
ref = np.isnan(linear_interpolation)
nearest_interpolation = anisotropy_function_nearest_interpolation(
points[ref,:])
linear_interpolation[ref] = nearest_interpolation
return linear_interpolation
r_test = np.linspace(0, 20, 100)
theta_test = np.linspace(0, 180, 100)
r_test_mesh, theta_test_mesh = np.meshgrid(
r_test, theta_test)
r_test_mesh_flat = r_test_mesh.ravel()
theta_test_mesh_flat = theta_test_mesh.ravel()
anisotropy_test = anisotropy_function(r_test_mesh_flat, theta_test_mesh_flat)
anisotropy_test = np.reshape(anisotropy_test, np.shape(r_test_mesh))
plt.contourf(r_test, theta_test, anisotropy_test, 1000)
c = plt.colorbar()
plt.xlabel("Radius (cm)")
plt.ylabel("Theta (degrees)")
c.set_label("Anisotropy function")
In [16]:
marker_cycle = itertools.cycle((',', 'o', '*', '^'))
colour_list = [
item['color']
for item in mpl.rcParams['axes.prop_cycle']
]
colour_cycle = itertools.cycle(colour_list)
theta_linspace = np.linspace(
np.array(np.min(anisotropy_function_data.theta_deg)),
np.array(np.max(anisotropy_function_data.theta_deg)), 1000)
for data in anisotropy_function_data.T:
colour = next(colour_cycle)
plt.scatter(
data.theta_deg, data, label="radius = {} cm".format(np.array(data.radius_cm)),
marker=next(marker_cycle), c=colour)
interpolation = anisotropy_function(np.array(data.radius_cm)*np.ones_like(theta_linspace), theta_linspace)
plt.plot(theta_linspace, interpolation, c=colour)
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Theta (degrees)')
plt.ylabel('Anisotropy function')
Out[16]:
In [17]:
r_test_mesh, theta_test_mesh = np.meshgrid(
anisotropy_radius_cm, anisotropy_theta_deg)
r_test_mesh_flat = r_test_mesh.ravel()
theta_test_mesh_flat = theta_test_mesh.ravel()
anisotropy_test = anisotropy_function(r_test_mesh_flat, theta_test_mesh_flat)
anisotropy_test = np.reshape(anisotropy_test, np.shape(r_test_mesh))
assert np.all(anisotropy_test - anisotropy_function_data_array == 0)
To better understand the geometry function see the following interactive demonstration:
The method here for calculating the geometry function is based on the following paper:
King, R. P., Anderson, R. S. and Mills, M. D. (2001), Geometry function of a linear brachytherapy source. Journal of Applied Clinical Medical Physics, 2: 69–72. doi:10.1120/jacmp.v2i2.2615
In [18]:
def calculate_beta(length, radius, theta):
# https://dx.doi.org/10.1120/jacmp.v2i2.2615
theta2 = np.arctan(
radius * np.sin(theta) /
(radius * np.cos(theta) - length/2))
AP = radius * np.sin(theta)
SA = radius * np.cos(theta) + length/2
SP = np.sqrt(AP**2 + SA**2)
beta = np.arcsin(
length * np.sin(theta2) /
SP)
return beta
def geometry_function(length, radius, theta):
theta = theta / 180 * np.pi
result = np.nan * np.ones_like(radius)
ref0 = (theta == 0) | (theta == np.pi)
ref1 = (theta != 0) & (theta != np.pi)
result[ref0] = (
1 / (radius[ref0]**2 - length**2 / 4))
beta = calculate_beta(
length, radius[ref1], theta[ref1])
result[ref1] = (
beta / (length * radius[ref1] * np.sin(theta[ref1])))
return np.abs(result)
# def geometry_function(length, radius, theta):
# return radius**-2
def normalised_geometry_function(length, radius, theta):
geometry = geometry_function(length, radius, theta)
ref = radius == 0
geometry_0 = geometry_function(length, np.array([0.000001]*np.sum(ref)), theta[ref])
geometry[ref] = geometry_0
return geometry / geometry_function(length, np.array([1]), np.array([90]))
In [19]:
def tg43(radius, theta):
initial_shape = np.shape(radius)
radius_flat = np.ravel(radius)
theta_flat = np.ravel(theta)
geometry = normalised_geometry_function(length, radius_flat, theta_flat)
radial = radial_function(radius_flat)
anisotropy = anisotropy_function(radius_flat, theta_flat)
result = dose_rate_constant * geometry * radial * anisotropy
result = np.reshape(result, initial_shape)
return result
In [20]:
def calc_on_grid(calc_grid_positions, dwell_positions, dwell_directions):
calc_grid_vector = calc_grid_positions[:,None,:] - dwell_positions[None,:,:]
radius = np.sqrt(
calc_grid_vector[:,:,0]**2 +
calc_grid_vector[:,:,1]**2 +
calc_grid_vector[:,:,2]**2)
calc_unit_vector = calc_grid_vector / radius[:,:,None]
dot_product = (
dwell_directions[None,:,0] * calc_unit_vector[:,:,0] +
dwell_directions[None,:,1] * calc_unit_vector[:,:,1] +
dwell_directions[None,:,2] * calc_unit_vector[:,:,2])
theta = np.arccos(dot_product) * 180 / np.pi
nan_ref = np.isnan(theta)
theta[nan_ref] = 90
tg43_dose = tg43(radius, theta)
return tg43_dose, radius
In [21]:
def determine_calc_grid_positions(calc_x, calc_y, calc_z):
calc_x_mesh, calc_y_mesh, calc_z_mesh = np.meshgrid(
calc_x, calc_y, calc_z)
calc_x_mesh_flat = calc_x_mesh.ravel()
calc_y_mesh_flat = calc_y_mesh.ravel()
calc_z_mesh_flat = calc_z_mesh.ravel()
calc_grid_positions = np.array([
[calc_x_mesh_flat[i], calc_y_mesh_flat[i], calc_z_mesh_flat[i]]
for i in range(len(calc_x_mesh_flat))
])
return calc_grid_positions, np.shape(calc_x_mesh)
In [22]:
def tg43_on_grid(calc_x, calc_y, calc_z,
reference_air_kerma_rate, dwell_times, dwell_positions, dwell_directions,
return_radius=False):
calc_grid_positions, shape = determine_calc_grid_positions(
calc_x, calc_y, calc_z)
tg43_dose, radius = calc_on_grid(
calc_grid_positions, dwell_positions, dwell_directions)
total_dose = np.sum(
tg43_dose * reference_air_kerma_rate * dwell_times[None, :], axis=1)
# too_close = np.any(radius < length/2, axis=1)
# total_dose[too_close] = np.nan
total_dose = np.reshape(total_dose, shape)
if return_radius:
result = (total_dose, np.reshape(radius, shape))
else:
result = total_dose
return result
In [23]:
testing_calc_x = np.array(0)
testing_calc_y = np.array(QA_along_away_data.y_cm)
testing_calc_z = np.array(QA_along_away_data.z_cm)
test_reference_air_kerma_rate = 1
test_dwell_times = np.array([1])
testing_dwell_positions = np.array([[0, 0, 0]])
testing_dwell_directions = np.array([[0, 0, 1]])
testing_tg43_dose, testing_radius = tg43_on_grid(
testing_calc_x, testing_calc_y, testing_calc_z,
test_reference_air_kerma_rate, test_dwell_times,
testing_dwell_positions, testing_dwell_directions, return_radius=True)
In [24]:
test_data = np.array(QA_along_away_data.T)
test_data[test_data > 1000] = np.nan
mpl.rcParams['contour.negative_linestyle'] = 'solid'
dose_difference = 100 * (testing_tg43_dose[:,0,:] - test_data) / test_data
levels = np.arange(-0.8,0.8,0.1)
CS = plt.contour(testing_calc_z, testing_calc_y, dose_difference, levels, colors='k', linewidths=0.2)
plt.clabel(CS,
inline=1,
fmt='%1.1f',
fontsize=9)
plt.contourf(testing_calc_z, testing_calc_y, dose_difference, np.arange(-0.8,0.8,0.01), extend='both', cmap='Spectral_r')
c = plt.colorbar()
c.set_label('Percent relative dose difference')
plt.xlabel('z (cm)')
plt.ylabel('y (cm)')
Out[24]:
In [25]:
np.nanmax(np.abs(dose_difference))
Out[25]:
In [26]:
reference = np.squeeze(testing_radius > 0.6)
np.nanmax(np.abs(dose_difference[reference]))
Out[26]:
In [27]:
import dicom
from glob import glob
from scipy.interpolate import make_interp_spline
dose_filename = glob("private/RD1*")[0]
plan_filename = glob("private/RP1*")[0]
# dose_filename = "test_dicom/RD_y5cm.dcm"
# plan_filename = "test_dicom/RP_y5cm.dcm"
In [28]:
dcm_dose = dicom.read_file(dose_filename, force=True)
dcm_plan = dicom.read_file(plan_filename, force=True)
In [29]:
dcm_plan.SourceSequence[0]
Out[29]:
In [30]:
reference_air_kerma_rate = np.float(dcm_plan.SourceSequence[0].ReferenceAirKermaRate) / 360000
reference_air_kerma_rate
Out[30]:
In [31]:
dcm_plan.ApplicationSetupSequence[0].ChannelSequence[0]
Out[31]:
In [32]:
def pull_dwells(dcm):
dwell_channels = []
all_dwell_positions = []
dwell_times = []
number_of_channels = len(dcm.ApplicationSetupSequence[0].ChannelSequence)
for i in range(number_of_channels):
ChannelSequence = dcm.ApplicationSetupSequence[0].ChannelSequence[i]
if len(ChannelSequence.dir("FinalCumulativeTimeWeight")) != 0:
BrachyControlPointSequence = (
ChannelSequence.BrachyControlPointSequence)
number_of_dwells_in_channel = len(BrachyControlPointSequence)//2
channel_time_weight = (
ChannelSequence.ChannelTotalTime / ChannelSequence.FinalCumulativeTimeWeight)
channel_dwells_positions = []
for j in range(0,number_of_dwells_in_channel*2,2):
if len(BrachyControlPointSequence[j].dir("ControlPoint3DPosition")) != 0:
channel_dwells_positions.append(
BrachyControlPointSequence[j].ControlPoint3DPosition)
dwell_channels.append(i+1)
assert (
BrachyControlPointSequence[j].ControlPoint3DPosition ==
BrachyControlPointSequence[j+1].ControlPoint3DPosition
)
uncorrected_dwell_time = (
BrachyControlPointSequence[j+1].CumulativeTimeWeight -
BrachyControlPointSequence[j].CumulativeTimeWeight
)
dwell_times.append(
uncorrected_dwell_time * channel_time_weight)
all_dwell_positions += channel_dwells_positions
dwell_channels = np.array(dwell_channels)
dwell_positions = np.array(all_dwell_positions).astype(float) / 10 # convert from mm to cm
dwell_times = np.array(dwell_times)
return dwell_positions, dwell_channels, dwell_times
dwell_positions, dwell_channels, dwell_times = pull_dwells(dcm_plan)
In [33]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
scatt = ax.scatter(
dwell_positions[:,0],
dwell_positions[:,1],
dwell_positions[:,2], c=dwell_times, s=dwell_times*15)
c = plt.colorbar(scatt)
c.set_label('Dwell time (s)')
ax.set(xlabel='x (cm)', ylabel='y (cm)', zlabel='z (cm)')
# plt.figure()
# plt.scatter(dwell_positions[:,0], dwell_positions[:,1])
Out[33]:
In [34]:
def pull_catheter_coords(dcm, channel):
index = channel - 1
# This line isn't correct for every dicom file. Specifically the "index" as given here
# doesn't uniquely define that given catheter channel. There is a specific tag within
# the contour sequence which should be used instead.
# The risk is that there is a single contour sequence that doesn't refer to a catheter
# instead I suspect that is where info such as catheter storage direction and the like
# is stored. It just so happens in this file that this private sequence is at the end
# in this example file. That is not always the case.
catheter_coords = dcm[0x300f,0x1000][0].ROIContourSequence[index].ContourSequence[0].ContourData
# mm to cm
x = np.array(catheter_coords[0::3]) / 10
y = np.array(catheter_coords[1::3]) / 10
z = np.array(catheter_coords[2::3]) / 10
return x, y, z
catheter_x, catheter_y, catheter_z = pull_catheter_coords(dcm_plan, 1)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
catheter_x,
catheter_y,
catheter_z)
ax.set(xlabel='x', ylabel='y', zlabel='z')
Out[34]:
In [35]:
print(dcm_plan[0x300f,0x1000][0].ROIContourSequence[0])
In [36]:
def dwell_direction_from_catheter_coords(catheter_x, catheter_y, catheter_z,
dwell_x, dwell_y, dwell_z):
t = np.linspace(0, 1, len(catheter_x))
spl = make_interp_spline(t, np.c_[
catheter_x, catheter_y, catheter_z
], k=1)
t_find_closest = np.linspace(0,1,10000)
x_find_closest, y_find_closest, z_find_closest = spl(t_find_closest).T
distance = np.sqrt(
(dwell_x[None,:] - x_find_closest[:,None])**2 +
(dwell_y[None,:] - y_find_closest[:,None])**2 +
(dwell_z[None,:] - z_find_closest[:,None])**2
)
dwell_params = t_find_closest[np.argmin(distance, axis=0)]
dwell_derivs = spl.derivative()(dwell_params).T
# dwell_directions = -dwell_derivs / np.linalg.norm(dwell_derivs, axis=0)
dwell_directions = dwell_derivs / np.linalg.norm(dwell_derivs, axis=0)
return dwell_directions.T
def determine_dwell_directions(dcm, dwell_positions, dwell_channels):
catheter_numbers = np.unique(dwell_channels)
dwell_directions = np.nan * np.ones_like(dwell_positions)
for catheter_number in catheter_numbers:
catheter_x, catheter_y, catheter_z = pull_catheter_coords(dcm, catheter_number)
current_ref = dwell_channels==catheter_number
current_dwell_positions = dwell_positions[current_ref,:]
current_dwell_directions = dwell_direction_from_catheter_coords(
catheter_x, catheter_y, catheter_z,
current_dwell_positions[:,0], current_dwell_positions[:,1],
current_dwell_positions[:,2])
dwell_directions[current_ref, :] = current_dwell_directions
return dwell_directions
dwell_directions = determine_dwell_directions(dcm_plan, dwell_positions, dwell_channels)
# Need to make a dicom file specific correction based off end/catheter or catheter/end
# private dicom tag
# dwell_directions = -dwell_directions
plt.plot(dwell_directions)
Out[36]:
In [37]:
def show_dwell_directions(catheter_x, catheter_y, catheter_z,
dwell_x, dwell_y, dwell_z, dwell_directions):
fig = plt.figure(figsize=(10, 10))
gs = gridspec.GridSpec(2, 2)
ax3 = fig.add_subplot(gs[1,0])
ax2 = fig.add_subplot(gs[0,1], sharex=ax3)
ax1 = fig.add_subplot(gs[0,0], sharey=ax2, sharex=ax3)
# ax1 = fig.add_subplot(gs[0,0])
# ax2 = fig.add_subplot(gs[0,1], sharey=ax1)
# ax3 = fig.add_subplot(gs[1,0], sharex=ax1)
ax4 = fig.add_subplot(gs[1,1], projection='3d')
ax1.plot(catheter_x, catheter_z, '--', alpha=0.2)
ax2.plot(catheter_y, catheter_z, '--', alpha=0.2)
ax3.plot(catheter_x, catheter_y, '--', alpha=0.2)
ax1.scatter(dwell_x, dwell_z)
ax2.scatter(dwell_y, dwell_z)
ax3.scatter(dwell_x, dwell_y)
ax1.quiver(
dwell_x,
dwell_z,
dwell_directions[:,0], dwell_directions[:,2],
pivot='mid', scale=10
)
ax2.quiver(
dwell_y,
dwell_z,
dwell_directions[:,1], dwell_directions[:,2],
pivot='mid', scale=10
)
ax3.quiver(
dwell_x,
dwell_y,
dwell_directions[:,0], dwell_directions[:,1],
pivot='mid', scale=10
)
ax1.set(aspect=1, xlabel='x (cm)', ylabel='z (cm)')
ax2.set(aspect=1, xlabel='y (cm)', ylabel='z (cm)')
ax3.set(aspect=1, xlabel='x (cm)', ylabel='y (cm)')
return ax4
test_catheter_x = np.arange(0, 3, 0.1)
test_catheter_y = test_catheter_x**2
test_catheter_z = 3*test_catheter_x
test_dwell_x = np.array([0.5, 1, 2, 2.5])
test_dwell_y = test_dwell_x**2
test_dwell_z = 3*test_dwell_x
test_dwell_directions = dwell_direction_from_catheter_coords(
test_catheter_x, test_catheter_y, test_catheter_z,
test_dwell_x, test_dwell_y, test_dwell_z)
show_dwell_directions(
test_catheter_x, test_catheter_y, test_catheter_z,
test_dwell_x, test_dwell_y, test_dwell_z, test_dwell_directions)
Out[37]:
In [38]:
def test_catheter(catheter_number):
catheter_x, catheter_y, catheter_z = pull_catheter_coords(dcm_plan, catheter_number)
min_z = np.min(dwell_positions[:,2]) - 1
max_z = np.max(dwell_positions[:,2]) + 1
catheter_ref = (catheter_z >= min_z) & (catheter_z <= max_z)
ax = show_dwell_directions(
catheter_x[catheter_ref], catheter_y[catheter_ref], catheter_z[catheter_ref],
dwell_positions[dwell_channels == catheter_number][:,0],
dwell_positions[dwell_channels == catheter_number][:,1],
dwell_positions[dwell_channels == catheter_number][:,2],
dwell_directions[dwell_channels==catheter_number,:])
ax.scatter(
dwell_positions[dwell_channels == catheter_number][:,0],
dwell_positions[dwell_channels == catheter_number][:,1],
dwell_positions[dwell_channels == catheter_number][:,2])
ax.scatter(
dwell_positions[dwell_channels != catheter_number][:,0],
dwell_positions[dwell_channels != catheter_number][:,1],
dwell_positions[dwell_channels != catheter_number][:,2])
ax.plot(
catheter_x[catheter_ref],
catheter_y[catheter_ref],
catheter_z[catheter_ref])
ax.set(xlabel='x', ylabel='y', zlabel='z')
catheter_numbers = np.unique(dwell_channels)
for catheter_number in catheter_numbers:
display(Markdown("#### Catheter Number {}".format(catheter_number)))
test_catheter(catheter_number)
plt.show()
In [39]:
# The x, y, and z defined here have not been sufficiently verified
# They do not necessarily match either what is within Dicom nor what is within your
# TPS. Please verify these and see if they are what you expect them to be.
# If these functions are incorrect or there is a better choice of dimension definitions
# please contact me by creating an issue within the github repository:
# https://github.com/SimonBiggs/npgamma/issues
# If you are able to validate these functions please contact me in the same way.
def load_dose_from_dicom(dcm):
"""Imports the dose in matplotlib format, with the following index mapping:
i = y
j = x
k = z
Therefore when using this function to have the coords match the same order,
ie. coords_reference = (y, x, z)
"""
pixels = np.transpose(
dcm.pixel_array, (1, 2, 0))
dose = pixels * dcm.DoseGridScaling
return dose
def load_xyz_from_dicom(dcm):
"""Although this coordinate pull from Dicom works in the scenarios tested
this is not an official x, y, z pull. It needs further confirmation.
"""
resolution = np.array(
dcm.PixelSpacing).astype(float)
# Does the first index match x?
# Haven't tested with differing grid sizes in x and y directions.
dx = resolution[0]
# The use of dcm.Columns here is under question
x = (
dcm.ImagePositionPatient[0] +
np.arange(0, dcm.Columns * dx, dx))
# Does the second index match y?
# Haven't tested with differing grid sizes in x and y directions.
dy = resolution[1]
# The use of dcm.Rows here is under question
y = (
dcm.ImagePositionPatient[1] +
np.arange(0, dcm.Rows * dy, dy))
# Is this correct?
z = (
np.array(dcm.GridFrameOffsetVector) +
dcm.ImagePositionPatient[2])
return x, y, z
In [40]:
tps_dose = load_dose_from_dicom(dcm_dose)
max_dose = np.max(tps_dose)
# Oncentra
estimated_perscription = max_dose / 10
relevant_dose_deviation = estimated_perscription * 0.03
max_dose_per_slice = np.max(np.max(tps_dose, axis=0), axis=0)
relevant_slice = max_dose_per_slice > estimated_perscription
In [41]:
x_raw, y_raw, z_raw = load_xyz_from_dicom(dcm_dose)
calc_grid_slice_skip = 5
calc_x = x_raw / 10
calc_y = y_raw / 10
calc_z = z_raw[relevant_slice] / 10
calc_z = calc_z[::calc_grid_slice_skip]
calc_grid_positions = determine_calc_grid_positions(calc_x, calc_y, calc_z)
In [42]:
tg43_dose = tg43_on_grid(
calc_x, calc_y, calc_z, reference_air_kerma_rate,
dwell_times, dwell_positions, dwell_directions)
In [43]:
tg43_dose[tg43_dose>max_dose] = max_dose
In [44]:
# Could be axis definition issues here
for z in calc_z:
tps_z_ref = np.where(z_raw/10 == z)[0][0]
calc_z_ref = np.where(calc_z == z)[0][0]
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(16,4))
c1 = ax1.contourf(
calc_x, calc_y, tg43_dose[:,:,calc_z_ref], np.arange(0,max_dose/2,0.1), extend='max')
ax1.set_title("Calculated dose @ z = {0:.2f}cm".format(z))
ax1.set_xlabel("x (cm)")
ax1.set_ylabel("y (cm)")
cbar = plt.colorbar(c1, ax=ax1)
cbar.ax.set_ylabel('Dose (Gy)')
c2 = ax2.contourf(
calc_x, calc_y, tps_dose[:,:,tps_z_ref], np.arange(0,max_dose/2,0.01), extend='max')
ax2.set_title("TPS dose @ z = {0:.2f}cm".format(z))
ax2.set_xlabel("x (cm)")
ax2.set_ylabel("y (cm)")
cbar = plt.colorbar(c2, ax=ax2)
cbar.ax.set_ylabel('Dose (Gy)')
dose_diff = tg43_dose[:,:,calc_z_ref] - tps_dose[:,:,tps_z_ref]
c3 = ax3.contourf(
calc_x, calc_y, dose_diff, np.arange(
-relevant_dose_deviation,relevant_dose_deviation,0.001), extend='both', cmap='Spectral_r')
ax3.set_title("Dose difference @ z = {0:.2f}cm".format(z))
ax3.set_xlabel("x (cm)")
ax3.set_ylabel("y (cm)")
cbar = plt.colorbar(c3, ax=ax3)
cbar.ax.set_ylabel('Dose Difference (Gy)')
plt.show()
In [45]:
# z = calc_z[len(calc_z)//2]
# tps_z_ref = np.where(z_raw/10 == z)[0][0]
# calc_z_ref = np.where(calc_z == z)[0][0]
# fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharey=True, figsize=(16,4))
# c1 = ax1.contourf(
# calc_x, calc_y, tg43_dose[:,:,calc_z_ref], np.arange(0,max_dose/2,0.001), extend='max')
# ax1.set_title("Calculated dose @ z = {0:.2f}cm".format(z))
# ax1.set_xlabel("x (cm)")
# ax1.set_ylabel("y (cm)")
# cbar = plt.colorbar(c1, ax=ax1)
# cbar.ax.set_ylabel('Dose (Gy)')
# dose_diff = tg43_dose[:,:,calc_z_ref] - tps_dose[:,:,tps_z_ref]
# c2 = ax2.contourf(
# calc_x, calc_y, dose_diff, np.arange(-relevant_dose_deviation,relevant_dose_deviation,0.0001), extend='both')
# ax2.set_title("Dose difference @ z = {0:.2f}cm".format(z))
# ax2.set_xlabel("x (cm)")
# ax2.set_ylabel("y (cm)")
# cbar = plt.colorbar(c2, ax=ax2)
# cbar.ax.set_ylabel('Dose Difference (Gy)')
# rel_dose_diff = 100 * dose_diff / tg43_dose[:,:,calc_z_ref]
# c3 = ax3.contourf(
# calc_x, calc_y, rel_dose_diff, np.arange(-3,3,0.0001), extend='both')
# ax3.set_title("Percent relative dose difference @ z = {0:.2f}cm".format(z))
# ax3.set_xlabel("x (cm)")
# ax3.set_ylabel("y (cm)")
# cbar = plt.colorbar(c3, ax=ax3)
# cbar.ax.set_ylabel('Percent Dose Difference')
# plt.show()